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We present a method for rapid modeling of new Bragg ultra-cold atom interferometer ( AI) designs 
useful for assessing the performance of such interferometers. The method simulates the overall effect 
on the condensate wave function in a given AI design using two separate elements. These are (1) 
modeling the effect of a Bragg pulse on the wave function and (2) approximating the evolution of 
the wave function during the intervals between the pulses. The actual sequence of these pulses and 
intervals is then followed to determine the approximate final wave function from which the interfer- 
ence pattern can be calculated. The exact evolution between pulses is assumed to be governed by 
the Gross-Pitaevskii (GP) equation whose solution is approximated using a Lagrangian Variational 
Method to facilitate rapid estimation of performance. The method presented here is an extension 
of an earlier one that was used to analyze the results of an experiment [J.E. Simsarian, et al., Phys. 
Rev. Lett. 83, 2040 (2000)], where the phase of a Bose-Einstein condensate was measured using 
a Mach-Zehnder-type Bragg AI. We have developed both ID and 3D versions of this method and 
we have determined their validity by comparing their predicted interference patterns with those 
obtained by numerical integration of the ID GP equation and with the results of the above experi- 
ment. We find excellent agreement between the ID interference patterns predicted by this method 
and those found by the GP equation. We show that we can reproduce all of the results of that 
experiment without recourse to an ad hoc velocity-kick correction needed by the earlier method, 
including some experimental results that the earlier model did not predict. We also found that this 
method provides estimates of ID interference patterns at least four orders-of-magnitude faster than 
direct numerical solution of the ID GP equation. 

PACS numbers: 03.75.Dg,67.85.Hj,03.67.Lx,03.75.Kk,42.50.Gy 
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I. INTRODUCTION 

It is possible to use ideas inspired by advances in 
quantum information science (QIS) to devise improved- 
performance matter-wave interferometers. A recent ex- 
ample of this has been seen in neutron interferometry 
where the idea of decoherence-free subspaces was used 
to redesign a neutron interferometer to reduce the effect 
of mechanical shaking on the interference contrast [1, 2]. 
The use of ideas from QIS to drive new neutron inter- 
ferometer designs may also be possible for atom interfer- 
ometers. Two promising areas of QIS where this could 
happen include decoherence avoidance and minimization. 
Interferometer applications where QlS-inspired redesigns 
may result in improved performance include precision 
navigation and metrology. This paper presents a tool 
for rapid assessment of new atom-interferometer designs 
for such applications. 

Atom interferometers (AI), where laser light is applied 
to ultra-cold atoms, have many applications. These in- 
clude quantum decoherence [3, 4], properties of Bose- 
Einstein condensates [5-8] , precision measurement of the 
fine-structure constant [9, 10] and testing the charge neu- 
trality of atoms [11]. Atom interferometers are also used 
in many precision measurement devices. These include 
gravimeters, gyroscopes, and gradiometers which all have 
important applications in precision navigation [12-14]. 
Atom interferometers also have applications in atomic 
physics such as atomic polarizability measurements and 



Casimir-Polder potentials for atoms near surfaces [15]. 
More uses of atom interferometry are described in Ref. 
[16]. 

With the advent of gaseous Bose-Einstein condensates 
(BEC) [17-21], strong interest has developed in using AIs 
for precision metrology [10, 22-28]. Most of these ultra- 
cold atom interferometers were of the standard Mach- 
Zehnder design. However, some more recent precision 
interferometers [10] have different designs. This also sug- 
gests that advances in interferometer design may lead to 
significant AI performance gains. 

There are many factors that can limit the performance 
of an AI. Some of these include mirror vibration, random 
initial motion of BECs at birth, stray light, external mag- 
netic fields, and errors in the frequency or intensity of the 
applied laser light. Atom interferometers confined on an 
atom chip can have other problems related to atom loss, 
heating, and decoherence [29] . One of the ways in which 
some of these factors may be addressed is with new AI 
designs. 

In order to pursue the program of drawing ideas from 
QIS to inspire new AI designs, it will be necessary to de- 
velop tools that can be used to provide rapid assessment 
of the performance of these new designs. In this work we 
present a method for rapid simulation of the behavior of 
condensates in Bragg interferometers. We assume that 
the evolution of the condensate between Bragg pulses is 
described by the Gross-Pitaevskii (GP) equation. The 
method approximates the evolution of the condensate 



wave function by modeling the effect of individual pulses 
and its evolution between pulses. Thus the final conden- 
sate wave function can be found enabling the prediction 
of the final interference pattern. As will be seen below, 
our method provides reasonable estimates of AI behav- 
ior in a time that is four orders-of~magnitude faster than 
that needed for numerical solution of the (ID) GP equa- 
tion to simulate a standard Mach-Zehndcr AI. The time- 
savings factor for 3D simulation will be greater. Such a 
tool will be essential for preliminary testing of new AI 
designs having more pulses and longer evolution times. 

In Section II we describe the standard 7r/2-7r-7r/2 
Bragg interferometer followed in Section III by a gen- 
eral overview of the two elements of our Bragg proto- 
typer model and the details of the original model used 
in the analysis of a Bragg interferometer experiment 
carried out at NIST [30]. In the original model each 
definite-momentum cloud in the condensate wave func- 
tion was represented by a single gaussian. An extra ve- 
locity "kick" , caused by the repulsion of the separating 
clouds after a 7r/2 Bragg pulse, was needed to obtain 
agreement with the experimental results. A study of this 
kick for a single pulse was conducted as a function of 
interaction strength of the condensate and the results 
are presented in Section IV. Section V presents the two- 
cloud version of our Bragg prototypcr model for both the 
ID and 3D cases. Finally, Section VI contains a summary 
and discussion of the possible applications of the method. 
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FIG. 1: (color online) This figure shows the sequence of Bragg 
pulses (represented by green arrows) applied to a condensate 
in a 7r/2-7r-7r/2 Bragg interferometer. The vertical (horizon- 
tal) direction represents space (time). See text for further 
details. 



cloud and starts the slow cloud so that the two clouds 
come back together. The final pulse is applied after a 
time T2 when the two clouds again overlap causing each 
cloud to split again. After a further time T3 there is a pair 
of overlapping slow clouds and a pair of overlapping fast 
clouds that give rise to the final pattern that is imaged 
at that moment. 



II. BRAGG ATOM INTERFEROMETER 

A Bose-Einstein condensate (BEG) can be coher- 
ently split into two clouds, a fast-moving cloud and a 
slow- moving cloud, through the application of a Bragg 
pulse [31]. If the condensate is stationary when a Bragg 
pulse is applied to it, then the result will be two clouds, 
one that remains stationary and another whose momen- 
tum is /lAk = 7i(ki — k2) where ki (k2) is the photon 
momentum of the higher (lower) frequency laser beam. 
If the condensate is moving with this momentum when 
the Bragg pulse is applied, then the result is again two 
clouds one of which keeps its momentum while the other 
cloud's momentum is reduced by ?iAk. In either case, 
the net result of applying a Bragg pulse to a condensate 
is a fast cloud and a slow cloud. One of these clouds has 
the momentum of the original cloud and the momentum 
of the other cloud is increased (decreased) if the original 
cloud was slow (fast). 

A Mach-Zehndcr-type Bragg interferometer can be 
constructed by applying three Bragg pulses in the se- 
quence 7r/2-7r-7r/2 with variable time intervals between 
them [32]. This is shown in Fig. 1 where the open circle 
shows the initial condensate which may be released from 
the trap and allowed to expand for a time Tg. The first 
7r/2 pulse splits the condensate into a slow cloud (upper 
path) and a fast cloud (lower path). After a time inter- 
val Ti the TT pulse is applied which stops the fast (lower) 



III. BRAGG INTERFEROMETER 
PROTOTYPER MODELS 

The operation of an arbitrary Bragg AI can be specified 
by stating the times and angles of the Bragg pulses ap- 
plied to the condensate. Although it is possible to apply 
a Bragg pulse having an arbitrary angle, 9, here we will 
restrict our attention to pulses where 9 = n ot 9 = 7r/2 
radians. Hence these pulses will either split clouds into 
two equal pieces, one fast and one slow, or will leave the 
condensate whole and merely swap its fast or slow veloc- 
ity. Thus we only consider interferometer sequences that 
are composed of tt/2 and/or tt pulses. 



A. Overview of prototyper models 

Our general prototyper model will enable us to ap- 
proximate the evolution of the condensate wave function 
through all of the steps of any given interferometer se- 
quence. We assume that the duration of all Bragg pulses 
is short compared to the characteristic time for collective 
effects of the condensate to be manifested. This time is 
usually h/ n where /z is the chemical potential of the con- 
densate. This means that we assume that any changes 
in the momentum space distribution of the condensate 
atoms caused by a Bragg pulse happen instantaneously. 
We further assume that the characteristic size of any mo- 



iTLcntuni change is large in the sense that the wavelength 
of the photon causing the change is small compared to 
the size of the condensate. Finally we assume that the 
evolution of the condensate between pulses is governed 
by the Gross-Pitaevskii equation. 

Our Bragg prototyper method therefore has two es- 
sential elements: (1) an approximation of the effect of 
a Bragg pulse on the condensate wave function, and (2) 
a model for approximating the GP-governed condensate 
wave function behavior between pulses. The method con- 
sists of applying these two elements to the pulse sequence 
of the given interferometer design to produce an approx- 
imate final condensate wave function so that predictions 
about the measured interference patterns can be made. 

For the first element, as will be described in more detail 
below, we represent the condensate wave function at any 
moment as a superposition of a number of fast or slow 
gaussian clouds. We model the effect of Bragg pulses as 
follows. If the wave function for a given cloud in the con- 
densate wave function is V'(r,^) before the Bragg pulse, 
then the change in this wave function is, for 
slow clouds: 

'i{jir,t) ^ ^(V'(r,t)-e-'''^e'^'^-""V'(r,t)) (7r/2 pulse) 
v2 



The GP equation is given by 



V'(r,i) ^ -e-"^e*'^''-Xr,i) (tt pulse) 



and for fast clouds: 



^(r,i) ^ -^(e**V(r,t) + e-*^''--^(r,i)) (^2 pulse) 
V'(r,t) ^ e*'^e-*'^*'-Xr, i) (tt pulse) (2) 

where Ak is the momentum change for Bragg pulses 
defined above. The factor </> is the phase of the mov- 
ing standing wave in the center of the initial atomic 
wavepacket in the middle of the Bragg pulse [32]. Thus 
the action of a 7r/2 pulse is to double the number of 
clouds, adding a new fast cloud on top of a previously 
existing slow cloud and adding a new slow cloud on top 
of a previously existing fast cloud. In each case the shape 
of the previously existing cloud is unchanged. The action 
of a TT pulse is to convert a previously existing slow (fast) 
cloud into a fast (slow) cloud. 

For the second element we use the Lagrangian Vari- 
ational Method (LVM) [33, 34] to approximate the con- 
densate evolution between pulses. Although it is possible 
to solve the 3D GP equation to determine this evolution, 
this is not practical for rapid estimation of the final con- 
densate wave function. The LVM provides approximate 
solutions to the GP equation in the form of equations 
of motion for time-dependent parameters that appear in 
an assumed trial wave function. Thus the exact solution 
of the GP equation that requires the solution of a 3-1-1 
partial differential equation is traded for approximate so- 
lutions that can be obtained by solving a system of or- 
dinary differential equations in time. We briefly review 
this method now. 
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a* 



-V2* + 14,,p(r)* + .giV|*|-* (3) 



dt 2m 

where we assume that any trapping potential is har- 
monic: 

Vtrap(r) == \mLolx^ + imw^y^ + Imuolz^ . (4) 

where m is the mass of a condensate atom. 

Here we will introduce scaled variables that will be 
used throughout the rest of the paper. First we choose a 
length unit appropriate to the harmonic potential: Lg = 
(?i/2TOCj)^/^ where Hi = [wxijJyiOzY^^ and then define the 
energy unit as Eq = h /2mL^ and the time unit as 
To = h/Eg. We then introduce scaled position and time 
variables as x = x/Lo,y = y/Lo,z = z/Lo,i = t/T^, 
and use barred quantities in general to represent quan- 
tities expressed in the scaled units. If, additionally, we 
write the condensate wave function in scaled units as 
^I* = ^/Lq'" then the GP equation becomes 



.a* 



-V^^ + Vtre.p{Y)^ + gNm * 



(5) 



(6) 



where V^ = d'^/dx'^ + d'^/df + d^/dz"^ and 

T> l-\ 1 2 -2 I 1 2 -2 I 1 2-2 

Vtrap(r) = jl,xX + ^-lyV + ^-1,Z 

and where ^ri = i^-n/w, rj = x,y,z and g = g/{EoLQ). 

The LVM produces equations of motion for the m 
time-dependent variational parameters that appear in a 
given trial wave function '!/'tiiai(f ; qi{t), . . . , qm{t))- The 
equations of motion are obtained from the LVM La- 
grangian via the usual Euler-Lagrange equations of mo- 
tion: 

The LVM Lagrangian, in turn, is computed by integrat- 
ing the LVM Lagrangian density 



^LVM 



few, 



(i))= /dV£[V;tdal(f,i)]. (8) 



Finally, the LVM Lagrangian density that corresponds to 
the GP equation is given by 

C[^] = i(V;^:-^>0+VV^*-W + Virap(f)|^|' 

+ \9N\^P\\ (9) 

where ij^t denotes the partial derivative of ip with respect 
to t. In order to get equations of motion relevant for a 
Bragg interferometer, we must choose a trial wave func- 
tion. In this work we present equations of motion for 
single-cloud gaussian trial wave functions in three di- 
mensions and two-cloud gaussian trial wave functions for 
both one and three dimensions. 

An alternative approach to efficient approximate solu- 
tion of the GP equation in interferometric applications 
has recently been presented by Jamison et al. [10]. It is 
based on a generalization of the Thomas-Fermi method 
originally proposed by Castin and Dum. [35] 
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FIG. 2: (color online) This figure shows the evolution of a ID 
condensate that is initially released and allowed to expand 
until T = 1 scaled time unit at which time a 7r/2 Bragg pulse 
is applied. The scaled velocity imparted by the Bragg pulse 
to the fast cloud is vl = 16 scaled velocity units, (a) The 
initial condensate density just after the trap is turned off, 
(b) after the tt/2 Bragg pulse is applied, the density exhibits 
rapid oscillations in the region where the fast and slow clouds 
overlap, (c) and (d) nearly complete and complete separation. 
The value of gN — 10 in scaled units. 



B. 3D single-cloud LVM model 

In the single-cloud LVM model we choose the trial 
wave function to have the form of a single, three- 
dimensional gaussian wavepacket, as was done previously 
in Rcf. [33] 



*(f,i)=^(i) n 



,-{n-no{i)f/2wl{i)+ia^{i]n+iPvit)n'' 



ri—x,y,z 



(10) 

Here (xo^yo, zq) are the coordinates of the center of 
the wavepacket. The quantities {ax,a.y,az) are the lin- 
ear phase coefficients which govern the motion of the 
wavepacket center. The {iVx,Wy,Wz) are the widths of 
the gaussian along the three axes and the [Px.Py, Bz) are 
the quadratic phase coefficients and govern the evolution 
of the widths. Finally, A{t) is a normalization coefficient 
that will be removed from the Lagrangian later when the 
normalization constraint is imposed. We can take A to 
be real because, if A had a phase, it would represent an 
overall wave function phase which would not be physical. 
The equations of motion are derived as described above 
and the result is [33] 
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^fjo - 2/3,,fjo T] = x,y,z. 



(11) 



It is worth noting that our equations of motion differ 



slightly from those in Ref. [33] because of differences in 
the definition of scaled units. 



NIST experiment and the single— cloud LVM 
model 



A Bragg AI was implemented at NIST and used to im- 
age the phase evolution of an evolving Bose-Einstein con- 
densate [30]. In that experiment, a condensate of about 
1.8 X 10^ sodium atoms was held in a magnetic trap with 
trapping frequencies oJx ~ V2ujy = IuJz = 27r x 27 Hz and 
subjected to -k JI-'k-'k jl pulse sequence. This sequence 
was performed both with the trap left on and with the 
trap turned off and the condensate allowed to expand for 
a time Tq. With the trap on, the time conditions were 
fixed so that Ti = T2 = T where T was typically 1-2 ms 
and the clouds were allowed to expand in order for them 
to separate before imaging. In the trap-off case, a series 
of runs was carried out in which Ti was held fixed at 1 
ms while T2 was varied such that the cloud overlap at the 
time of the final Bragg pulse ranged from fast cloud just 
arriving at the slow cloud until it had passed through and 
was just leaving. This series of runs was performed for 
an expansion time of Tq = 1 ms and repeated for Tq — A 
ms. 

The results of this experiment were analyzed using the 
single-cloud LVM just described [30]. While the results 
of the single-cloud LVM model agreed well with exper- 
iment for runs performed with the trap on, it did not 
agree with experiment for trap-off cases. This discrep- 
ancy was due to the presence of an extra relative veloc- 
ity between interfering clouds at the moment of the final 
Bragg pulse. The extra velocity was caused by repulsion 
between overlapping clouds which occurred just after the 
first 7r/2 pulse and again just before the second 7r/2 pulse. 
Agreement between theory and experiment was achieved 
by adding in by hand a small relative velocity correc- 
tion to the condensate wave function predicted by the 
single-cloud model. It is clear that any model able to 
account for this repulsion would need to include at least 
two clouds. In order to derive a simple model we need to 
validate a key approximation by studying this correction 
for a single 7r/2 Bragg pulse. 



IV. REPULSION STUDY FOR SINGLE n/2 
BRAGG PULSE 



Before turning to a two-cloud model, we studied the 
effects of fast/slow cloud repulsion on the final relative 
velocity of the separating clouds after a single 7r/2 Bragg 
pulse. We performed this study by simulating the ap- 
plication of such a pulse on a ID condensate and its 
subsequent evolution by numerical solution of the GP 
equation. 

In what follows we shall give the value of quantities in 
ID scaled units which is a special case of the 3D scaled 



units given above. In these units, all q uantities are de- 
fined in terms of the length unit Lp = ^J'h/2raujQ which, 
in turn, is tied to a reference frequency, wq. This fre- 
quency is often the trap frequency but need not be as 
in the case where the trap has been turned off. Since 
these quantities are scaled out of the problem, it will be 
useful to give a numerical example of the sizes of the 
scaled units. Thus, given a quasi-lD ^^Rb condensate 
confined in a a; = 27r x 10 Hz trap, the length unit is 
Lq = 2.4 /Ltm, the time unit is Tq = 15.8 ms, and the ve- 
locity unit Vq = 0.015 cm/s. The value of the interaction 
strength is varied over the range < gmN < 200 so that 
the transition from non-interacting up to the Thomas- 
Fermi regime could be studied. Here, N is the number 
of condensate atoms. 

Figure 2 shows a typical simulation where a condensate 
is released from the trap (panel (a)) and allowed to ex- 
pand for until T — 1 scaled time unit at which time a 7r/2 
Bragg pulse is applied splitting the condensate into fast 
(on the right) and slow clouds. During the separation the 
two clouds push each other apart so that the fast cloud 
moves with a velocity Vf = vl + Sv that is slightly larger 
than the recoil velocity vl caused by the laser light and 
the slow cloud drifts backwards with velocity —Sv. In SI 
units, the recoil velocity is vl ~ hAk/m and, in scaled 
units it can be expressed as vl = vl/{Lq/Tq) = 2/Sk. 
The interaction strength of the initial condensate in this 
example was gN = 10 in scaled units. 

To study this velocity "kick", 5v, we used the GP sim- 
ulations to determine its value as a function of the in- 
teraction strength of the initial condensate. The value 
of Sv for a given value of gN was obtained by running 
a simulation where a 7r/2 Bragg pulse was applied. The 
velocities of the fast and slow clouds were determined, for 
each value of gN, as follows. In each run the two clouds 
were allow to separate fully after the pulse; a midpoint 
between the two clouds, Xmid, and a time after which 
both clouds were fully separated, isep, were then deter- 
mined; and then the expectation value of x was computed 
numerically for each cloud separately at each time step 
in the range i > tgcp '■ 
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(12) 



where L is the length of the numerical grid used in the 
GP simulation. Care was taken to make sure that none 
of the clouds got close to the edges of the grid. Finally, 
the velocities of the fast and slow clouds were extracted 
by fitting straight lines to the Xsiow and Xfast results to 
obtain Vsiow and Ufast- The extra velocity "kicks" for 
each cloud due to repulsion were determined by (5wfast = 
^fast — Vl and (Jugiow = Vsiow Convergence runs of the 
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FIG. 3: (color online) This figure shows the velocity "kick", 
5v, versus interaction strength, gN . The open circles are 
the results obtained from GP-equation simulations; the three 
curves show estimates of the velocity kick obtained by equat- 
ing the difi'erence in total kinetic energy of the interacting 
system and the non-interacting system with 2/3 of the total 
interaction energy. These estimates are derived from (1) the 
exact GP initial state (solid curve), (2) LVM-approximate ini- 
tial state (dashed curve) and (3) Thomas-Fermi approximate 
initial state (dotted curve). 



GP solver for finer space and time steps showed that 
(^"Wsiow = Sviast = 5v over the entire interaction strength 
range. This shows that momentum was conserved and 
supports the picture of two equal-mass clouds pushing 
against each other as they separate. 

Fig. 3 contains a graph of the velocity kick versus in- 
teraction strength for < gN < 200. For this set of 
simulations, the Bragg pulse was applied and the trap 
was turned off simultaneously at t ~ 0. For the results 
shown, the recoil velocity was wl = 10 scaled velocity 
units. In the example mentioned above this would give 
^^Rb atoms a recoil velocity of 0.15 cm/s. The result 
displayed in Fig. 3 can be quantitatively understood in a 
simple way as we now explain. 

The velocity kick, 6v, is the result of mutual repulsion 
between the fast and slow clouds as they separate after 
the Bragg pulse. From a classical viewpoint, this repul- 
sion will result in a change in the total kinetic energy 
of the center-of-mass (CM) motion of the two clouds. 
Thus the difference between the total CM kinetic energy 
of the two interacting clouds and the kinetic energy of the 
two non-interacting clouds should be equal to the energy 
available for repulsion assuming the repulsion produces 
no distortion. 

Since the Bragg pulse splits the A^-atom condensate 
into two equal pieces we can express this equality as (we 
begin in SI units and then convert): 



c/„ 



^) ^m{vL + SvY + {^) ^miSv) 



, 2 J 



\mvl 



(13) 



where C/rcp is the many-body energy available for differ- 
ent clouds to repel each other. We can derive an expres- 
sion for this from the total many-body interaction energy 
given by [20] 



Uint^hNiN-l)g dx\^{x,t)\ 



(14) 



where *I'(a;,i) is the condensate wave function after the 
7r/2 Bragg pulse and has the general form 



^(a;,i) =Mx,t) 



lAfea 



1p2ix,t). 



(15) 



We have chosen to distinguish between tpi and V'2, even 
though they are the same in our model, for bookkeeping 
purposes that will become apparent below. Substituting 
Eq. (15) into (14) we have 



Un 



fiV(iV-l) 



dx 



i^ii*+i^2i^+4ivinv'2 



+ 2 (iViP + IV2P) ^i^se^^'^^ + (^i>2) 



^2 2iAkx 



(16) 



Here we will make a crucial approximation. This ap- 
proximation will be tested by comparison of the velocity 
kicks predicted here with those determined by numerical 
solution of the GP equation and will be used again in 
deriving a two-cloud LVM model. 

We assume here that Afc is large enough so that all of 
the integrals containing exponentials such as exp(±iAfca;) 
and exp(±2iAfca;) in Eq. (16) can be neglected. This is 
equivalent to assuming that the wavelength of the Bragg 
pulse light is small compared to the size of the conden- 
sate. Hence terms having these exponentials oscillate 
rapidly and the integrals containing them approximately 
average to zero. The result of this approximation is that 
all of the terms inside the curly braces in Eq. (16) can be 
neglected and we can write 

/CXD 
-00 
= C/sclf,l + t^sclf,2 + C/rcp. (17) 

This last expression suggests a picture of the evolution 
of fast and slow clouds during separation. This picture 
depends on two assumptions: (1) the separating clouds 
do not distort significantly so that the form of the wave 
function in Eq. (15) is maintained and, (2) the wavevec- 
tor, Ak, is large enough so that the approximation in Eq. 
(17) is valid. In this case, the fast/slow cloud evolution 
during separation divides into three categories: (1) self 
interaction of the fast cloud, (2) self interaction of the 
slow cloud, and (3) fast/slow cloud interaction. The en- 
ergy available for this last interaction is suggested by the 



above equation: 

/"OO 

C/,.ep = igiV(iV-l) / A\^i\^\ij2\^dx 



(18) 



igiV^ / l^fdx = |C/int. 



Here we assume that A^ ^ 1 and that ipi ^ ip2 ^ 
where ^ is the (unit norm) condensate wave function 
before application of the Bragg pulse as in Eqs. (1) and 

(2). 

This fact also leads to the last equality because, from 
Eq. (17), it is clear that the total energy of interaction 
is one part self interaction of cloud 1, one part self in- 
teraction of cloud 2, and 4 parts cloud-cloud interaction. 
The self interaction terms, according to our picture, are 
energies available for expansion while the cloud-cloud in- 
teraction either distorts the cloud shapes and/or changes 
the velocities of their centers-of-mass. We assume no 
distortion so all of this energy is assumed available for 
giving the clouds a velocity kick. We can now derive this 
kick by substituting the approximate expression for C/j-op 
into Eq. (13) and canceling common factors of N/2: 

gN ItPl'^dx = ^m{vL+Svf + ^m{5vf 



^mvl 



gN 



'dx = l(vL+Svy + \iSv) 



l„-i2 



(19) 



where in second line above we have converted back to 
scaled units. Thus we can now write an expression for 
the velocity kick: 



'5'" = [j^L + 2Usp] - 5WL, 



where 



Usp = gN / li/'l dx 



(20) 



(21) 



and vl — 2Ak. We can think of the separating clouds 
being pushed apart by a spring in between them and Ugp 
is the initial energy stored in the spring. 

The comparison of values of Sv predicted by the above 
model with those determined from numerical solution of 
the GP equation are shown in Fig. 3. The discrete points 
are the numerically determined values and the remaining 
three curves arc computed via Eq. (20) where Ugp has 
been calculated using three expressions for ip, the initial 
state condensate wave function. These three versions of 
■0 were (1) the exact initial state from the GP simula- 
tion (solid curve), (2) an LVM gaussian cloud where the 
gaussian width was the stationary value for a ID conden- 
sate confined in the initial harmonic trap (dashed line), 
and (3) the Thomas-Fermi-approximatc solution of the 
GP equation for the trapped condensate (dotted line). 



All three expressions, as can be seen from the graph, 
gave good estimates. We found that using the exact GP 
initial state gave near-perfect agreement (to four deci- 
mal places) with the kick determined from numerical GP. 
The LVM gaussian performed very well for small gN and 
less well for large while the Thomas-Fermi works well for 
large gN and less well for small. 

This agreement between the numerically determined 
kicks and the estimates from our heuristic model lends 
support to our picture of the effect of cloud/cloud inter- 
action during separation. We will use this in what follows 
to develop a two-cloud LVM technique for modeling a full 
Bragg interferometer. 



V. TWO-CLOUD LVM MODEL 

In this section we present a two-cloud LVM model to 
approximate the evolution of the condensate wave func- 
tion following a 7r/2 Bragg pulse. The presence of two 
clouds will enable the model to account for cloud-cloud 
interactions and should be able to predict the velocity 
kick accurately. Below we present both ID and 3D ver- 
sions of this model. This will enable us to make quanti- 
tative comparisons of our model with the results of ID 
GP simulations of the entire 7r/2-7r-7r/2 interferometer. 
The 3D model will be useful for assessment of real-world 
AI designs. In addition, we also present comparisons of 
our 3D two-cloud model with the results of the original 
NIST experiment. 



A. ID two-cloud LVM 

The trial wave function for the ID two-cloud LVM 
model is taken to be a sum of two ID gaussian wavepack- 
ets _where one of them is boosted to a velocity of vl = 
2Afc: 



f'(x,i) 



V2 



/g/lCS,*) _^ giAfcSg/2(2,t)A 



(22) 



where 



fjix,t) 



{x-Xj{t)Y 



+ iaj{t)x + i^{t)x'^ (23) 



where j = 1,2. We note here that both the slow cloud 
(cloud 1) and the fast cloud share the same width, w, and 
quadratic phase curvature, /3 but have differing centers 
and linear phase coefficients. This assumption is borne 
out in multiple GP simulations as can be seen, for ex- 
ample, in Fig. 2. Imposing the normalization condition 
yields the following constraint: \Aid\'^ivtt^^'^ = 1 where 
all terms containing exponentials such as e^*^*^^ were 
neglected. 

Carrying out the procedure described in Section III A, 
we calculate the LVM Lagrangian by inserting the trial 



wave function into the Lagrangian density and integrat- 
ing. If we neglect rapidly oscillating terms and impose 
the normalization constraint we obtain the following re- 
sult. 

Lid = \aixi A- \a2X2 + \P {x\ + x\ + up') 

v2 



1 

2w^ 



i(ai 



2/3Si 



5 ("2 



2/3x2 



2„T.2 



+ (Afc) (52 + 2px2) + i (Afc) + 2p''w 



1~,2 



+ ^7 SI+S^+i*^ 



+ \A(2l)^'^w] \ 



1-K2e" 



-(xi— a;2)"/2to" 



(24) 



where 7 ~ i^/ojq is the ratio of the actual trap frequency 
to the frequency used to define the length unit. 

We find the equations of motion using the ordinary 
Euler-Lagrange equations (see Eq. (7)). After some re- 
arrangement these equations can be expressed as follows: 



Xi 



■7^x1 



i^i^o^ (xi,X2,U') 



a^2 + 7 X2 



12 
"■^12 



(X1,X2,W) 
i^i,l^)(xi,X2,«)) 



(25) 



where 

F12 [XI,X2,W) 



and 

Fi,i^)(xi,X2,u;) 



__2gN__ 
(27r)i/2iB3 



(xi 



X2)e-(^i-*2''/2'^' 
(26) 



-( m. "i 

V(27r)i/2iDV 



-{X1-X2) /2w^ 



(27) 

where ¥[2 can be thought of as the "force" of repulsion 
between the separating clouds and Fw roughly thought 
of as the "force" causing the width to change. 

The above equations constitute a closed system of 
equations to be solved for Xi, xi, X2 , xi, w, and w. 
Once these equations are solved, the remaining parame- 
ters appearing in the trial wave function can be obtained 
as follows. 







ai = \xi - 


-2Pxi 


a.2 = \X2 - 


-2^X2 



Afc 



(28) 



The above system of equations have some interesting 
properties that are analogous to those for Newton's sec- 
ond law. For example, if we introduce the center of mass 
of the two-cloud system, Xcm = (Si -|-X2)/2, and the rel- 
ative coordinate, Xrei = xi — X2 we find their equations 
of motion by adding and subtracting the first two of Eqs. 
(25) respectively. They are 



Xcm I 7 Xcm 
Xrel ~r 7 Xrel 





2F^l '{xi,X2,w) 



2F, 



{ID) 
12 
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(29) 
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FIG. 4: (color online) These plots show a comparison of 
the interference patterns resulting from a 7r/2-7r-7r/2 Bragg 
AT computed by the GP equation (dashed line) and by the 
Bragg prototyper method. The interferometer times were 
(in scaled time units) To = i, Ti — 4, T2 variable, and 
Ts = 13.4. The times for T2 are (clockwise from upper 
left) T2 = 3.6, 3.8, 4.25, 4.5 scaled time units. The interac- 
tion strength for all runs was gN — 10 scaled units. 



where, in the last equality, we noted that F^fj o^^Y 
depends on the relative coordinate and the width. The 
first of the above equations can be solved immediately by 
inspection and is equivalent to the classical result that the 
CM motion of a system only depends on external forces. 
This reduces the total number of equations that must be 
solved numerically to the equation for w in Eqs. (25) and 
the equation for Xrei above. 

There is also a conserved "energy" that can be written 
as 



EiD 



\xl + \x\ + jW^ + Uid{xItX2tW) 



(30) 



where the "potential energy" Uid is given by 

1 2 -2 , 1 2 -2 , 1 2-2 



Tt I- - -\ 12-2,12-2,12-2, 2 

UiD{Xi,X2,w) = 57 Xi 4- 57 X2 + 57 W +^ 



(31) 

Equations (25) can all be written as the second time 
derivative of each coordinate equals the negative partial 
derivative of the above potential energy with respect to 
the corresponding coordinate (e.g., xi = —dUiD/dxi). 
This constant of the motion can be used to estimate the 
velocity kick and is also useful as a check on numerics. 

Equations (25) and (28) can be used with the rest of 
the Bragg prototyper model to simulate Bragg AI be- 
havior. We have used this ID Bragg prototyper model 
to predict the results of a 7r/2-7r-7r/2 Bragg AI and have 
also simulated such an AI with the ID GP equation. In 
these runs, the condensate was released from the trap and 



allowed to expand for Tq = 1 scaled time unit. The 7r/2 
Bragg pulse was then applied and, after a time interval 
T2 = 4 scaled units, a tt pulse was applied followed by a 
variable time interval, T2, at which time the final Bragg 
pulse was turned on. The clouds were allowed to evolve 
for an additional Tt, ~ 13.4 time units. Figure 4 shows 
the comparison of the interference patterns predicted by 
the Bragg prototyper model and by the GP equation. 
Each graph corresponds to a different value of T2. In the 
figure, the values were (clockwise from upper left panel) 
T2 = 3.6,3.8,4.25,4.5 time units. 

The GP equation was solved using a Crank-Nicolson 
algorithm on a space grid of width L = 800 length units 
and divided into N^ = 32768 space steps over a total time 
of imax = 22 scaled time units using Nt = 1200000 time 
steps. The initial state was computed by integrating the 
GP equation in imaginary time with the trap on. The 
time needed to obtain fully converged final interference 
patterns was about five hours of run time on a commodity 
laptop. The time required to obtain the model results on 
the same computer was about 1 second, approximately 
18,000 times faster. We expect that the speedup factor 
for 3D solutions to be 2-3 orders-of-magnitude greater. 
This makes our model an essential tool for assessment 
of more complicated atom interferometers which result 
from new AI designs inspired by quantum information 
science. 

It is easy to see that the Bragg prototyper model re- 
produces the number and spacing of the GP-equation 
interference fringes in all cases. The major difference is 
that the width and height of the GP fringes are wider and 
shorter, respectively, than those in the LVM pattern. It 
is likely that the area under the curve in corresponding 
GP and LVM fringes is the same. This would make wider 
fringes lower and would imply that the number of atoms 
in each GP fringe was approximately equal to the number 
of atoms in the corresponding LVM fringe. 



B. 3D two-cloud LVM 

We also derived a 3D version of the two-cloud LVM 
model. The derivation is a straightforward generalization 
of the ID version with all of the same assumptions and 
approximations. We include the highlights of its deriva- 
tion for completeness. 

The trial wavefunction for the 3D two-cloud LVM 
model is a sum of two 3D gaussians: 



*(f,i) 



A. 
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(32) 



where 



F,(f,t)= E {-^^S^+'i^^vV + ^f)] (33) 



2zi;2 



and J = 1,2. In three dimensions we allow for the pos- 
sibility of differing widths, iD^, 77 = x,y,z along the 
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FIG. 5: (color online) This figure exhibits a comparison of the results of the 3D two-cloud Bragg AI model with the experimental 
results given in Fig. 2 of Ref. [30] where the NIST experiment discussed earlier in the text was described. In that experiment, 
the evolving phase of a Bose-Einstein condensate wave function was probed using a 7r/2-7r-7r/2 Bragg interferometer, (a) 
Interference patterns obtained for interferometer runs (top row experiment, bottom row theory) where To = 4 ms, Ti — 1 
ms, Ts = 2 ms, and T2 was varied so that the cloud-center spacings (Sx) at the last Bragg pulse were the values given in the 
figure; (b) A quantitative comparison of the case for 5x — llfim (second from right end in panel (a)). The theory curve was 
normalized to match the highest experimental peak. 



three axes, however we assume that these widths are the 
same for both clouds. This holds for the correspond- 
ing phase curvature coefficients Pn, "H = x,y,z as well. 
It is straightforward to calculate the normalization con- 
straint as \A3£)\'^WxWyWyTr^^^ = 1. We note again that 
it was necessary to neglect rapidly oscillating terms that 
contained exponentials such as e^'^*^'' to arrive at this 
result. 

The Lagrangian for this trial wave function thus be- 
comes 



L^n — 



Y^ I \oiin-ni + \a2nm + ^A, (?7? + ?72 + w^) 

+ 2^7 + 2^r>' + hi {fi +nl + wl) 



U(2< 



) ' WxWyWs; 
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T] = x,y,z 



2wl 



(34) 



where 7^ = LiJr,/u) and w is the geometric average of the 
three trap frequencies that is used in the definition of the 
3D length unit. 

The equations of motion that arise from the above La- 
grangian are 



W' 



m + lr,Vl = 


F^S^^v^^-v^,^) 


m + i^m = 


-F^^'f (fi,f2,w) 


\ + ll^r^ = 


i^(;^,,^)(fi,f2,w) 



r\ = x,y,z. (35) 



where r^ = (xj , y^ , Zj ) with j = 1,2 and w = {wx ,Wy,Wz) 
and where the "force" terms on the right-hand-sides 



above are given by 
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(36) 



The rest of the equations connect the widths, and CM 
positions and velocities to the parameters that actually 
appear in the trial wave function. These are 



Oi2ri 



4wn 
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(37) 



The 3D version of the 2-cloud model shares some of the 
properties associated with the ID version. These include 
(1) the motion of the CM, fcm ~ (f 1 + F2)/2 is governed 
only by the trapping force and, (2) there is a conserved 
energy given by 



E^D = \ Yl {fi+vl+wl]+UM'ri,r2,^) (38) 
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FIG. 6: (color online) (a) The spatial fringe frequency versus the cloud separation at the instant the final 7r/2 pulse is applied. 
The data is taken from Fig. 3 of Ref. [30]. The steeper line corresponds to To = 4 ms, and the other line to To = 1 ms. (b) The 
relative velocity between interfering clouds versus T, the interferometer time. Data is taken from Fig. 5 of the cited paper, (c) 
An interference pattern comparison between 3D one-cloud LVM (bottom), experiment (middle), and 3D two-cloud LVM (top), 
with the trap on during all the steps of the interferometer and T = 1 ms. 



where 
C73D(fi,f2,w) = Y, {^7^('7? + ^^+wf, 
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These equations can now be used as a part of the two- 
gaussian-cloud LVM method to predict the behavior of 
atom interferometers. We have applied this model to the 
NIST experiment described earlier and detailed in Ref. 
[30]. Figure 5(a) shows interference patterns obtained 
in the experiment (top row) compared with the model. 
The timings of interferometer runs that produced these 
patterns were Tg = 4 ms, Ti ~ 1 ms, and T3 = 2 ms. The 
values of T2 were varied so that the spacing of the cloud 
centers, at the moment the final Bragg pulse is applied, 
6x, was as shown in Fig. 5(a). 

We see that there is good agreement with the ex- 
perimental patterns in terms of number and spacing of 
fringes. Figure 5(b) shows a more quantitative compari- 
son for the case of Sx = ll/itm. The data was taken from 
Fig. 2 of Ref. [30] and the theory is the result of the 3D 
two-cloud Bragg AI prototyper. All of these experiments 
were carried out with the trapping potential turned off. 
Since the kick is included naturally in our model, there 
was no need to add a velocity kick correction to achieve 
this level of agreement as there was for the one-cloud 
model. 

In the NIST experiment, the spatial fringe frequency, 
K, at the last Bragg pulse was measured as a function 
of 5x for the trap off case. The comparison of these re- 
sults with the 3D two-cloud LVM is shown in Fig. 6(a) 



where the data is taken from Fig. 3 of Ref. [30]. In 
our LVM method, the atom density at the time of the 
last pulse oscillates as cos((a22; — c>:ix)x). Thus we have 
klvm = Oi2x — o.ix evaluated at the time of the last pulse. 
The above comparison shows excellent agreement with 
experiment. The relative velocity between the clouds in 
a given pair was also measured for different values of 
Ti = T2 = T when the trap was left on. The data from 
Fig. 5(b) of Ref. [30] is shown in Fig. 6(b) along with 
the results of the LVM for the relative velocity versus T. 
Again we find good agreement with the experiment. 

The Fig. 6(c) shows a comparison of the 3D one-cloud 
LVM (bottom), 3D two-cloud LVM (top), and experi- 
ment (middle), for the interference pattern resulting from 
an interferometer run in which the trap was on and T ~ 1 
ms. Although the one-cloud and two-cloud LVM inter- 
ference patterns are both qualitatively similar to the ex- 
perimental one, it is clear that the two-cloud pattern 
agrees better. Thus the apparent agreement between the 
one-cloud LVM and experiment presented in Ref. [30] in 
the trap-on case without the correction was fortuitous. 



VI. DISCUSSION 

In this paper we have presented a method suitable for 
rapid estimation of interference patterns deriving from 
ultra-cold Bragg atom interferometers. The method 
achieves this by representing the condensate wave func- 
tion as a superposition of fast and slow gaussian clouds 
and then modeling the changes caused by (1) Bragg n/2 
and TT pulses and (2) by approximating the GP-equation 
evolution of the wave function during the intervals be- 
tween the pulses. Thus, by following the sequence of 
Bragg pulses and intervals that occur in a particular 
Bragg AI, it is possible to approximate rapidly the final 
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condensate wave function and thus calculate the expected 
interference pattern. 

In this model, Bragg pulses are assumed to change 
the wave function by instantaneous shifts in momentum 
space while rapid estimation of the effect of GP evolu- 
tion is approximated using an LVM technique. We have 
validated the ID version of this method by comparing its 
results with ID numerical simulations of Mach-Zehnder- 
type atom interferometers using the GP equation. The 
3D version was validated by comparison with experimen- 
tal results presented in Ref. [30]. We found that the 
method provides good agreement with interference pat- 
terns in both the ID and 3D cases as regards the number 
and spacing of fringes. In the ID case, we found that our 
model obtained the final interference patterns at least 
10,000 times faster than direct GP-equation simulation 
on the same computer. We expect this factor to be sev- 
eral orders-of-magnitude greater for the 3D case. 

While the comparison of the results of the method with 
ID GP simulations and with 3D experiment is quite good, 
there are some differences. In the ID case, we found 
that it was necessary to shift our interference patterns 
over to line up fringe locations with the GP results. The 
shifts were small; they were at most 10 scaled length units 
which was a small difference compared to the width of a 
fringe. We did not have enough information about the 
experiment to determine if this was necessary in the 3D 
case. Also, as can be seen in Fig. 4(a), the fringe heights 
differ from the GP pattern. This difference is also present 
in the 3D comparison. 

There are several possible reasons for these differences. 
First, our model tacitly assumes that separating clouds 
are not distorted by their mutual interaction since we 
model them as gaussians. In the ID case, this assumption 
holds reasonably but not perfectly well over the range of 
gN considered. GP simulations with larger and larger 
gN values show that clouds become more and more dis- 
torted as the interaction strength increases. Thus some 
of the energy available for repulsion can be diverted from 
changing the CM velocity into cloud distortion. 

Another reason for the difference may lie in our treat- 
ment of interacting clouds between the final Bragg pulse 
(which creates two fast and two slow clouds) and the time 
that an image is taken. During this interval, we use the 
two-cloud model to propagate each fast/slow pair sep- 
arately and combine them at the end. This neglected 
the interaction of the overlapping clouds (the two fast 
clouds overlap and the two slow clouds overlap) as the 
two cloud pairs move apart. In fact, the T3 interval is 
usually longer than all of the other intervals combined so 
there is ample opportunity for the two clouds in a pair 
to repel each other. This could be accounted for in a 
four-cloud model. 

We also presented a study of the extra velocity kick 
that the slow and fast clouds receive when a 7r/2 Bragg 
pulse is applied to a condensate. When the condensate 
splits, the two clouds push against each other as they 
separate. This causes the fast cloud to acquire an extra 



velocity, 5v, in addition to the velocity, v^, imparted by 
the light. The slow cloud recoils at velocity Sv due to 
conservation of momentum. This study was conducted 
in ID where the evolution of the condensate was sim- 
ulated with the GP equation and the value of 5v was 
determined directly from the simulation as a function of 
the interaction strength gN. The values of gN ranged 
from non-interacting up to well into the Thomas-Fermi 
regime. 

We found that these velocity kicks could be precisely 
predicted by setting the energy available for repulsion 
equal to the total kinetic energy of the interacting system 
minus the kinetic energy of the non-interacting system. 
We approximated the energy available for repulsion from 
the expression for the total interaction energy by neglect- 
ing rapidly oscillating terms involving the wavevector of 
the Bragg pulse light. The success of this procedure in re- 
producing the velocity kicks reinforces the picture of the 
interaction energy being partitioned into self-interaction 
energy of individual clouds and energy for cloud-cloud 
interaction. And furthermore, that the cloud-cloud in- 
teraction only produces a CM velocity change. 

It is also important to mention that there are some 
phenomena that occur in Bragg atom interferometers 
that the model described in this paper might be mod- 
ified to handle. Bragg processes often result in elastic 
scattering into intially unoccupied transverse momentum 
modes [36]. While these processes are not treated by 
the GP equation, it is possible to modify the effect of a 
pulse in these cases by adding clouds that occupy these 
momentum modes and neglecting their interaction with 
the mother condensate during separation. These modi- 
fications may also address other wave mixing processes 
that sometimes occur. Stray light can also cause prob- 
lems for Bragg atom interferometers especially for con- 
densates trapped near an atom chip. The major effect of 
stray light in these cases is to change the internal energy 
state of the atom. The model described above can be 
generalized to account for multiple internal levels of the 
atoms. 

Another difficulty for Bragg interferometers is imper- 
fections in the Bragg pulse wavefront which causes the 
laser intensity to vary across the condensate. In this case, 
different condensate atoms experience different pulse ar- 
eas which causes the "angle" of the Bragg pulse to vary. 
That is, not all atoms will experience either a tt or tt/2 
pulse. Our model can also be generalized to handle Bragg 
pulses of arbitrary angle. These "angle" errors in the 
pulses are analogous to imperfections in the thickness of 
the blades in a neutron interferometer. Such errors can 
be addressed by a new AI design that implements this 
QIS idea of the "power of one qubit" [37]. The imple- 
mentation of this idea for Bragg interferometers will be 
the subject of a forthcoming article. 

Atom interferometry with Bose-Einstein condensates 
hold the promise for applications in ultra-sensitive navi- 
gation and precision metrology. We noted earlier the idea 
that, because of the connection between quantum algo- 
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rithms and multi-particle intcrfcroiTLctcrs, there is great 
potential for using the advances in QIS to inspire ad- 
vances in precision interferometer designs. The Bragg AI 
prototyping method presented here represents a new tool 
for the rapid assessment of new Bragg AI designs. In the 
future we intend to apply our tool to design AIs that use 
QIS ideas for dccohcrcncc avoidance (e.g., decohercnce- 
frce subspaccs, [38]) and minimization (e.g., the power of 
one qubit,[37]). 
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